# -*- coding: utf-8 -*-
"""
Creating a normalised model from the original model

@author: obiri
"""
import mpctools as mpc
import numpy as np
import casadi as cs
import time

##16 states, 21 parameters. 
#nb: some of the remaining parameters are fixed measurements, one is state dependent

lb = 0.5 # lower bound of normalized model: lb*xs
ub = 1.5 # upper bound of normalized model: ub*xs
width = ub-lb

# xs = np.array([3.27367434e+10, 8.28831113e+10, 2.71741538e+03, 5.78604334e+00,
#         6.37656737e+01, 4.83137928e+00, 6.23866477e+01, 3.00039503e+03,
#         2.61960613e+09, 6.63233551e+09, 2.17448574e+03, 4.63001301e+00,
#         5.10255242e+01, 3.86608729e+00, 4.99220220e+01, 2.97694060e+03,
#         3.70349577e+01, 1.50000000e+01, 4.99220218e+01])

xs = np.array([3.27367434e+10, 8.28831113e+10, 2.71741538e+03, 5.78604334e+00,
        6.37656737e+01, 4.83137928e+00, 6.23866477e+01, 
        2.61960613e+09, 6.63233551e+09, 2.17448574e+03, 4.63001301e+00,
        5.10255242e+01, 3.86608729e+00, 4.99220220e+01, 
        3.70349577e+01, 4.99220218e+01])




num_fixed_param = 18
pn=np.zeros((num_fixed_param))
pn[0]= 1.76
pn[1]= 9.6e-03/60  #fix  #0.00016
pn[2]= 0.75
pn[3]= 0.075
pn[4]= 28.48
pn[5]= 171.76
pn[6]= 2.0 #fix
pn[7]= 0.45
pn[8]= 2
pn[9]= 2.6e8
pn[10]= 8e8
pn[11]= 4
pn[12]= 0.06/60      #0.001
pn[13]= 6.59e-10/60  #1.0983333333333334e-11
pn[14]= 5e5
pn[15]= 1560.0
pn[16]= 1.244
pn[17]= 4e2

# # ##Generating the true parameters for the plant by adding some bias
# np.random.seed(5)
# q=np.zeros((num_fixed_param))
# # q = np.random.normal(loc=q, scale=std_dev)    
# q[0] = np.random.normal(loc=pn[0], scale=0.0001)   
# q[1] = np.random.normal(loc=pn[1], scale=0.000001)   
# q[2] = np.random.normal(loc=pn[2], scale=0.00001) 
# q[3] = np.random.normal(loc=pn[3], scale=1e-6)   
# q[4] = np.random.normal(loc=pn[4], scale=1e-3) 
# q[5] = np.random.normal(loc=pn[5], scale=1e-3)   
# q[6] = np.random.normal(loc=pn[6], scale=1e-5) 
# q[7] = np.random.normal(loc=pn[7], scale=1e-5) 
# q[8] = np.random.normal(loc=pn[8], scale=1e-4) 
# q[9] = np.random.normal(loc=pn[9], scale=1e4) 
# q[10] = np.random.normal(loc=pn[10], scale=1e4) 
# q[11] = np.random.normal(loc=pn[11], scale=1e-3) 
# q[12] = np.random.normal(loc=pn[12], scale=1e-8) 
# q[13] = np.random.normal(loc=pn[13], scale=1e-11) 
# q[14] = np.random.normal(loc=pn[14], scale=1e2) 
# q[15] = np.random.normal(loc=pn[15], scale=1e-1) 
# q[16] = np.random.normal(loc=pn[16], scale=1e-3) 
# q[17] = np.random.normal(loc=pn[17], scale=1e-3) 
# ###################################################################


factor = 1.01*np.ones(18)
factor[0] = 1.1
factor[3] = 1.1
factor[12] = 1.1
factor[13] = 1.1
factor[15] = 1.1
q = factor*pn




##Plant Model uses nominal paramters + bias
def reactor_tank_uncert(x, u, w):   
    
    """
    An integrated model including upstream and buffer tank only.
    :param x: State vector (19,1)
    :param u: Input vector (8,1)
    :return: ODEs (19,1)
    """
    dxdt = cs.SX.zeros(16)
    
    # x = x_aug[0:16]

      
    #defining new variables for the normalised states
    x1 = x[0]*(width*xs[0])+lb*xs[0];      x2 = x[1]*(width*xs[1])+lb*xs[1];      x3 = x[2]*(width*xs[2])+lb*xs[2]; 
    x4 = x[3]*(width*xs[3])+lb*xs[3];      x5 = x[4]*(width*xs[4])+lb*xs[4];      x6 = x[5]*(width*xs[5])+lb*xs[5]; 
    x7 = x[6]*(width*xs[6])+lb*xs[6];      x8 = x[7]*(width*xs[7])+lb*xs[7];      x9 = x[8]*(width*xs[8])+lb*xs[8];
    x10 = x[9]*(width*xs[9])+lb*xs[9];     x11= x[10]*(width*xs[10])+lb*xs[10];   x12 = x[11]*(width*xs[11])+lb*xs[11];
    x13 = x[12]*(width*xs[12])+lb*xs[12];  x14= x[13]*(width*xs[13])+lb*xs[13];   x15 = x[14]*(width*xs[14])+lb*xs[14];
    x16 = x[15]*(width*xs[15])+lb*xs[15];  
    
    # p1 =  p[0]*(width*p_n[0])+lb*p_n[0];   p2 =  p[1]*(width*p_n[1])+lb*p_n[1];       p3 =  p[2]*(width*p_n[2])+lb*p_n[2]
    # p4 =  p[3]*(width*p_n[3])+lb*p_n[3];   p5 =  p[4]*(width*p_n[4])+lb*p_n[4]


    # Inputs of upstream model:
    F_in = u[0]    # inlet flow rate (L/min)
    F_1 = u[1]     # outlet flow rate of bioreactor (L/min)
    F_r = u[2]     # u[1]-u[3] # u[2]      # recycle flow rate (L/min)
    F_2 = u[3]     # Outlet flow rate of separator (L/min)
    GLC_in = u[4]  # inlet glucose concentration (mM)
    GLN_in = u[5]  # inlet glutamine concentration (mM)
    Tc = u[6]      # coolant temperature (d C)
    # Inputs of buffer tank:
    F_out_bf = u[7]  # Outlet flow rate of the buffer tank (L/min)  
    F_in_bf = u[3]   # Inlet flow rate of the buffer tank (L/min)
    #c_in_bf = x[14] # Inlet concentration of mAb (mg/L)
    c_in_bf = x14    ##changed to the scaled value of x[14]
    
    
    K_damm = q[0]#1.76     # ammonia constant for cell death (mM)
    K_dgln = q[1]#9.6e-03/60  # constant for glutamine degradation (min^(-1))
    K_glc =  q[2]#0.75      # Monod constant for glucose (mM)
    K_gln =  q[3] #0.075     # Monod constant for glutamine (mM)
    KI_amm = q[4]#28.48    # Monod constant for ammonia (mM)
    KI_lac = q[5]#171.76   # Monod constant for lactate (mM)
    m_glc =  4.9e-14/60   ###--ALWAYS FIXED--##   maintenance coefficients of glucose (mmol/cell/min) 
    n =      q[6]#2             # constant represents the increase of specific death as ammonia concentration increases (-)
    Y_ammgln = q[7]#0.45   # yield of ammonia from glutamine (mmol/mmol)
    Y_lacglc = q[8]#2.0    # yield of lactate from glucose (mmol/mmol)
    Y_Xglc =   q[9]#2.6e8    # yield of cells on glucose (cell/mmol)
    Y_Xgln =   q[10]#8e8      # yield of cells on glutamine (cell/mmol)
    alpha1 =   3.4e-13/60   ###--ALWAYS FIXED--##   #  # constants of glutamine maintenance coefficient (mM L/cell/min)
    alpha2 =   q[11]#4        # constants of glutamine maintenance coefficient (mM)
    # mu_max = (0.0016*x[16] - 0.0308)*5/60 # 0.058 * 5  # maximum specific growth rate (min^(-1))
    mu_max = (0.0016*x15 - 0.0308)*5/60    ##--ALWAYS FIXED--##   # 0.058 * 5  ##x[16] changed to the scaled value of x16
    #mu_dmax = (-0.0045*x[16] + 0.1682)/60 # 0.06    # maximum specific death rate (min^(-1))
    mu_dmax = q[12]# 0.06/60
    m_mabx =  q[13]#6.59e-10/60 # constant production of mAb by viable cells (mg/Cell/min) = Qmab_max
    
    # Add new parameters -- Ben
    Delta_H = q[14]#5e5     # 4.4e4# 1e-8
    rho = q[15]# 1560.0      # density of mixture is assumed to be density of glucose [g/L]
    cp =  q[16]#1.244        # specific heat capacity of mixture is assumed to be that of glucose [J/g/dC]
    U =   q[17]#4e2
    T_in =  37.0    ###--ALWAYS FIXED--##
    
    ####---these parameters were excluded because they are very easy to know, they willnot change even after a very long time----####
    # Parameters of buffer tank
    D = 5  # Diameter of the tank (dm)
    L = D  # Diameter of the tank (dm)
    Ac= np.pi*(D/2)**2  # Cross-sectional area (dm^2)
    V = Ac*L  # Volume of the tank (L)
           
    # Assumptions: 92% cell recycle rate and 20% mab retention
    Xvr = 0.92 * x1 * F_1 / F_r  # concentration of viable cells in recycle stream(cell/L)
    Xtr = 0.92 * x2 * F_1 / F_r  # total concentration of cells in recycle stream(cell/L)
    GLCr = 0.2 * x3 * F_1 / F_r  # glucose concentration in recycle stream(mM)
    GLNr = 0.2 * x4 * F_1 / F_r  # glutamine concentration in recycle stream(mM)
    LACr = 0.2 * x5 * F_1 / F_r  # lactate concentration in recycle stream(mM)
    AMMr = 0.2 * x6 * F_1 / F_r  # ammonia concentration in recycle stream(mM)
    mAbr = 0.2 * x7 * F_1 / F_r  # mAb concentration in recycle stream (mg/L)
            
    # ODEs
    f_lim = (x3 / (K_glc + x3)) * (x4 / (K_gln + x4))
    f_inh = (KI_lac/ (KI_lac + x5)) * (KI_amm / (KI_amm + x6))
    mu = mu_max * f_lim * f_inh
    mu_d = mu_dmax / (1 + ( K_damm / x6) ** n)
    
    dxdt[0]= (mu * x1 - mu_d * x1 - (F_in / 3.00039503e+03) * x1 + (F_r / 3.00039503e+03) * (Xvr - x1) +w[0]) /(width*xs[0])
    dxdt[1] = (mu * x1 + (F_r / 3.00039503e+03) * (Xtr - x2) - (F_in / 3.00039503e+03) * x2 +w[1]) /(width*xs[1])
    
    Q_glc = mu / Y_Xglc + m_glc
    dxdt[2] = (-Q_glc * x1 + (F_in / 3.00039503e+03) * (GLC_in - x3) + (F_r / 3.00039503e+03) * (GLCr - x3) +w[2]) /(width*xs[2])
    
    m_gln = (alpha1 * x4) / (alpha2 + x4)
    Q_gln = mu / Y_Xgln + m_gln
    dxdt[3] = (-Q_gln * x1 - K_dgln * x4 + (F_in / 3.00039503e+03) * (GLN_in - x4) + (F_r / 3.00039503e+03) * (GLNr - x4) +w[3]) /(width*xs[3])
    
    Q_lac = Y_lacglc * Q_glc
    dxdt[4] = (Q_lac * x1 - (F_in / 3.00039503e+03) * x5 + (F_r / 3.00039503e+03) * (LACr - x5) +w[4]) /(width*xs[4])
    
    Q_amm = Y_ammgln * Q_gln
    dxdt[5] = (Q_amm * x1 + K_dgln * x4 - (F_in / 3.00039503e+03) * x6 + (F_r / 3.00039503e+03) * (AMMr - x6) +w[5]) /(width*xs[5]) 
    dxdt[6] = (m_mabx * x1 - (F_in / 3.00039503e+03) * x7 + (F_r / 3.00039503e+03) * (mAbr - x7) +w[6]) /(width*xs[6])
    
    # dxdt[7] = (F_in + F_r - F_1)/(width*xs[7])
    # dxdt[7] = 0
    
    dxdt[7] = ((F_1 / 2.97694060e+03) * (x1 - x8) - (F_r / 2.97694060e+03) * (Xvr - x8) +w[7]) /(width*xs[7])
    dxdt[8] = ((F_1 / 2.97694060e+03) * (x2 - x9) - (F_r / 2.97694060e+03) * (Xtr - x9) +w[8]) /(width*xs[8])
    dxdt[9] = ((F_1 / 2.97694060e+03) * (x3 - x10) - (F_r / 2.97694060e+03) * (GLCr - x10) +w[9]) /(width*xs[9])
    dxdt[10] = ((F_1 / 2.97694060e+03) * (x4 - x11) - (F_r / 2.97694060e+03) * (GLNr - x11) +w[10]) /(width*xs[10])
    dxdt[11] = ((F_1 / 2.97694060e+03) * (x5 - x12) - (F_r / 2.97694060e+03) * (LACr - x12) +w[11]) /(width*xs[11])
    dxdt[12] = ((F_1 / 2.97694060e+03) * (x6 - x13) - (F_r / 2.97694060e+03) * (AMMr - x13) +w[12]) /(width*xs[12])
    dxdt[13] = ((F_1 / 2.97694060e+03) * (x7 - x14) - (F_r / 2.97694060e+03) * (mAbr - x14) +w[13]) /(width*xs[13])
    # dxdt[15] = (F_1 - F_2 - F_r) /(width*xs[15])
    # dxdt[15] = 0
    
    # Temperature dynamics -- Ben
    # Average mass of an mAb cell is 150 kDa = 2.4908084e-19 g
    dxdt[14] = ((F_in / 3.00039503e+03) * (T_in - x15) + (Delta_H / (rho  * cp)) * mu * x1 * 2.4908084e-19 + (U / (3.00039503e+03 * rho * cp)) * (Tc - x15) +w[14]) /(width*xs[14]) 
        
    # ODEs of buffer tank:
    dhdt = 1/Ac*(F_in_bf-F_out_bf)
    dcdt = F_in_bf/(Ac*1.50000000e+01)*(c_in_bf-x16)
    # dxdt[17] = dhdt /(width*xs[17])
    # dxdt[17] = 0
    dxdt[15] = (dcdt +w[15]) /(width*xs[15])
                                                                                                                                                         
    #dxdt = [dxdt1, dxdt2, dxdt3, dxdt4, dxdt5, dxdt6, dxdt7, dxdt8, dxdt9, dxdt10, dxdt11, dxdt12, dxdt13, dxdt14, dxdt15, dxdt16, dxdt17, dxdt18]
    return dxdt #(returns scaled values since this is now the scaled 



###ESTIMATOR MODEL which uses the nominal parameters
def reactor_tank(x, u, w):   
    
    """
    An integrated model including upstream and buffer tank only.
    :param x: State vector (19,1)
    :param u: Input vector (8,1)
    :return: ODEs (19,1)
    """
    dxdt = cs.SX.zeros(16)
    
    #defining new variables for the normalised states
    x1 = x[0]*(width*xs[0])+lb*xs[0];      x2 = x[1]*(width*xs[1])+lb*xs[1];      x3 = x[2]*(width*xs[2])+lb*xs[2]; 
    x4 = x[3]*(width*xs[3])+lb*xs[3];      x5 = x[4]*(width*xs[4])+lb*xs[4];      x6 = x[5]*(width*xs[5])+lb*xs[5]; 
    x7 = x[6]*(width*xs[6])+lb*xs[6];      x8 = x[7]*(width*xs[7])+lb*xs[7];      x9 = x[8]*(width*xs[8])+lb*xs[8];
    x10 = x[9]*(width*xs[9])+lb*xs[9];     x11= x[10]*(width*xs[10])+lb*xs[10];   x12 = x[11]*(width*xs[11])+lb*xs[11];
    x13 = x[12]*(width*xs[12])+lb*xs[12];  x14= x[13]*(width*xs[13])+lb*xs[13];   x15 = x[14]*(width*xs[14])+lb*xs[14];
    x16 = x[15]*(width*xs[15])+lb*xs[15];  
    
    
    
    #inputs remained the same as defined above
    # x8 = xs[7]; 
    # x16 = xs[15];
    # x18 = xs[17];
    
    # Inputs of upstream model:
    F_in = u[0]    # inlet flow rate (L/min)
    F_1 = u[1]     # outlet flow rate of bioreactor (L/min)
    F_r = u[2]     # u[1]-u[3] # u[2]      # recycle flow rate (L/min)
    F_2 = u[3]     # Outlet flow rate of separator (L/min)
    GLC_in = u[4]  # inlet glucose concentration (mM)
    GLN_in = u[5]  # inlet glutamine concentration (mM)
    Tc = u[6]      # coolant temperature (d C)
    # Inputs of buffer tank:
    F_out_bf = u[7]  # Outlet flow rate of the buffer tank (L/min)  
    F_in_bf = u[3]   # Inlet flow rate of the buffer tank (L/min)
    #c_in_bf = x[14] # Inlet concentration of mAb (mg/L)
    c_in_bf = x14    ##changed to the scaled value of x[14]
    

    #parameters - one of the parameters is state dependent so the parameters had to be brought below the normalised variables
    K_damm = 1.76     # ammonia constant for cell death (mM)
    K_dgln = 9.6e-03/60  # constant for glutamine degradation (min^(-1))
    K_glc = 0.75      # Monod constant for glucose (mM)
    K_gln = 0.075     # Monod constant for glutamine (mM)
    KI_amm = 28.48    # Monod constant for ammonia (mM)
    KI_lac = 171.76   # Monod constant for lactate (mM)
    m_glc = 4.9e-14/60   # maintenance coefficients of glucose (mmol/cell/min)
    n = 2             # constant represents the increase of specific death as ammonia concentration increases (-)
    Y_ammgln = 0.45   # yield of ammonia from glutamine (mmol/mmol)
    Y_lacglc = 2.0    # yield of lactate from glucose (mmol/mmol)
    Y_Xglc = 2.6e8    # yield of cells on glucose (cell/mmol)
    Y_Xgln = 8e8      # yield of cells on glutamine (cell/mmol)
    alpha1 = 3.4e-13/60  # constants of glutamine maintenance coefficient (mM L/cell/min)
    alpha2 = 4        # constants of glutamine maintenance coefficient (mM)
    # mu_max = (0.0016*x[16] - 0.0308)*5/60 # 0.058 * 5  # maximum specific growth rate (min^(-1))
    mu_max = (0.0016*x15 - 0.0308)*5/60 # 0.058 * 5  ##x[16] changed to the scaled value of x16
    #mu_dmax = (-0.0045*x[16] + 0.1682)/60 # 0.06    # maximum specific death rate (min^(-1))
    mu_dmax =0.06/60
    m_mabx = 6.59e-10/60 # constant production of mAb by viable cells (mg/Cell/min) = Qmab_max
    
    # Add new parameters -- Ben
    Delta_H = 5e5     # 4.4e4# 1e-8
    rho = 1560.0      # density of mixture is assumed to be density of glucose [g/L]
    cp = 1.244        # specific heat capacity of mixture is assumed to be that of glucose [J/g/dC]
    U = 4e2
    T_in = 37.0
    
    # Parameters of buffer tank
    D = 5  # Diameter of the tank (dm)
    L = D  # Length of the tank (dm)
    Ac = np.pi*(D/2)**2  # Cross-sectional area (dm^2)
    V = Ac*L  # Volume of the tank (L)
    
    
    # Assumptions: 92% cell recycle rate and 20% mab retention
    Xvr = 0.92 * x1 * F_1 / F_r  # concentration of viable cells in recycle stream(cell/L)
    Xtr = 0.92 * x2 * F_1 / F_r  # total concentration of cells in recycle stream(cell/L)
    GLCr = 0.2 * x3 * F_1 / F_r  # glucose concentration in recycle stream(mM)
    GLNr = 0.2 * x4 * F_1 / F_r  # glutamine concentration in recycle stream(mM)
    LACr = 0.2 * x5 * F_1 / F_r  # lactate concentration in recycle stream(mM)
    AMMr = 0.2 * x6 * F_1 / F_r  # ammonia concentration in recycle stream(mM)
    mAbr = 0.2 * x7 * F_1 / F_r  # mAb concentration in recycle stream (mg/L)
    
    
    # ODEs
    f_lim = (x3 / (K_glc + x3)) * (x4 / (K_gln + x4))
    f_inh = (KI_lac / (KI_lac + x5)) * (KI_amm / (KI_amm + x6))
    mu = mu_max * f_lim * f_inh
    mu_d = mu_dmax / (1 + (K_damm / x6) ** n)
    
    dxdt[0]= (mu * x1 - mu_d * x1 - (F_in / 3.00039503e+03) * x1 + (F_r / 3.00039503e+03) * (Xvr - x1) +w[0]) /(width*xs[0])
    dxdt[1] = (mu * x1 + (F_r / 3.00039503e+03) * (Xtr - x2) - (F_in / 3.00039503e+03) * x2 +w[1]) /(width*xs[1])
    
    Q_glc = mu / Y_Xglc + m_glc
    dxdt[2] = (-Q_glc * x1 + (F_in / 3.00039503e+03) * (GLC_in - x3) + (F_r / 3.00039503e+03) * (GLCr - x3) +w[2]) /(width*xs[2])
    
    m_gln = (alpha1 * x4) / (alpha2 + x4)
    Q_gln = mu / Y_Xgln + m_gln
    dxdt[3] = (-Q_gln * x1 - K_dgln * x4 + (F_in / 3.00039503e+03) * (GLN_in - x4) + (F_r / 3.00039503e+03) * (GLNr - x4) +w[3]) /(width*xs[3])
    
    Q_lac = Y_lacglc * Q_glc
    dxdt[4] = (Q_lac * x1 - (F_in / 3.00039503e+03) * x5 + (F_r / 3.00039503e+03) * (LACr - x5) +w[4]) /(width*xs[4])
    
    Q_amm = Y_ammgln * Q_gln
    dxdt[5] = (Q_amm * x1 + K_dgln * x4 - (F_in / 3.00039503e+03) * x6 + (F_r / 3.00039503e+03) * (AMMr - x6) +w[5]) /(width*xs[5]) 
    dxdt[6] = (m_mabx * x1 - (F_in / 3.00039503e+03) * x7 + (F_r / 3.00039503e+03) * (mAbr - x7) +w[6]) /(width*xs[6])
    
    # dxdt[7] = (F_in + F_r - F_1)/(width*xs[7])
    # dxdt[7] = 0
    
    dxdt[7] = ((F_1 / 2.97694060e+03) * (x1 - x8) - (F_r / 2.97694060e+03) * (Xvr - x8) +w[7]) /(width*xs[7])
    dxdt[8] = ((F_1 / 2.97694060e+03) * (x2 - x9) - (F_r / 2.97694060e+03) * (Xtr - x9) +w[8]) /(width*xs[8])
    dxdt[9] = ((F_1 / 2.97694060e+03) * (x3 - x10) - (F_r / 2.97694060e+03) * (GLCr - x10) +w[9]) /(width*xs[9])
    dxdt[10] = ((F_1 / 2.97694060e+03) * (x4 - x11) - (F_r / 2.97694060e+03) * (GLNr - x11) +w[10]) /(width*xs[10])
    dxdt[11] = ((F_1 / 2.97694060e+03) * (x5 - x12) - (F_r / 2.97694060e+03) * (LACr - x12) +w[11]) /(width*xs[11])
    dxdt[12] = ((F_1 / 2.97694060e+03) * (x6 - x13) - (F_r / 2.97694060e+03) * (AMMr - x13) +w[12]) /(width*xs[12])
    dxdt[13] = ((F_1 / 2.97694060e+03) * (x7 - x14) - (F_r / 2.97694060e+03) * (mAbr - x14) +w[13]) /(width*xs[13])
    # dxdt[15] = (F_1 - F_2 - F_r) /(width*xs[15])
    # dxdt[15] = 0
    
    # Temperature dynamics -- Ben
    # Average mass of an mAb cell is 150 kDa = 2.4908084e-19 g
    dxdt[14] = ((F_in / 3.00039503e+03) * (T_in - x15) + (Delta_H / (rho * cp)) * mu * x1 * 2.4908084e-19 + (U / (3.00039503e+03 * rho * cp)) * (Tc - x15) +w[14]) /(width*xs[14]) 
    
    
    # ODEs of buffer tank:
    dhdt = 1/Ac*(F_in_bf-F_out_bf)
    dcdt = F_in_bf/(Ac*1.50000000e+01)*(c_in_bf-x16)
    # dxdt[17] = dhdt /(width*xs[17])
    # dxdt[17] = 0
    dxdt[15] = (dcdt +w[15]) /(width*xs[15])
                                                                                                                                                         
    #dxdt = [dxdt1, dxdt2, dxdt3, dxdt4, dxdt5, dxdt6, dxdt7, dxdt8, dxdt9, dxdt10, dxdt11, dxdt12, dxdt13, dxdt14, dxdt15, dxdt16, dxdt17, dxdt18]
    return dxdt #(returns scaled values since this is now the scaled version of the ode )


num_measurements=7
num_states=16
def measurement_xy_model(x):            
    C_matrix = np.zeros((num_measurements, num_states))
    C_matrix[0,7]=1
    C_matrix[1,8]=1
    C_matrix[2,9]=1
    C_matrix[3,10]=1
    C_matrix[4,11]=1
    C_matrix[5,12]=1
    C_matrix[6,15]=1
    
    y = np.dot(C_matrix, x)
    return y





# ###------FOR PARAMETER ERROR PLOT AND RMSE CALCULATION----####

Nsim = 50
num_fixed_param = 18
p_nom = np.zeros((Nsim, num_fixed_param))

for i in range(pn.shape[0]):
    p_nom[:,i] = pn[i]
    

p_true = np.zeros((Nsim, num_fixed_param))
for i in range(q.shape[0]):
    p_true[:,i] = q[i]



def calculate_param_unscaled_all(p_nom, p_true):
    error = ((p_nom - p_true) / p_true )**2 # Element-wise relative error
    rmse = np.sqrt(np.mean(error, axis=1))  # RMSE for each simulation
    nrmse = np.mean(rmse)  # Mean of RMSE across simulations
    return nrmse
nrmse_p1 = calculate_param_unscaled_all(p_nom, p_true)
print("NRMSE:", nrmse_p1)



# def rmse_param_normalized_all(x_act, x_pred):
#     error = ((x_pred - x_act))**2 # Element-wise relative error
#     rmse = np.sqrt(np.mean(error, axis=1))  # RMSE for each simulation
#     nrmse = np.mean(rmse)  # Mean of RMSE across simulations
#     return nrmse
# rmse_value_p_normalized =rmse_param_normalized_all()






#########-------Plot for parameter ERROR TRAJECTORY------########
##The error is calculated based on the normalised x and xhat values
#The error is summed across columns to get the error across the states for 1 time step(sum_err)
#The trajectory is the total error for each time step, so it's a plot of 50 values
from matplotlib import pyplot as plt
Delta = 60 # Time step
tplot = np.arange(Nsim+1)*Delta 
err_p = np.abs(p_nom - p_true)

fig_0, axs = plt.subplots(1, figsize=(10,10)) 
sum_err_p = err_p.sum(axis=1)
axs.plot(tplot[:-1],sum_err_p, '-.b' )
axs.set_xlim([0, 800])
plt.ylabel('Error between nominal and true parameters')
plt.xlabel('Time (min)')
# np.save('sum_Perror1' , sum_err_p)


sum_Perror1 = 0.45*np.ones(50)
np.save('sum_Perror1' , sum_Perror1)




